Enhanced Algorithm for Desktop PC-Based Micro-Finite Element Modeling of Whole-Section Stiffness from in Vivo MR Images

Introduction: High-resolution MR image-based micro finite-element modeling (μFE) is a valuable tool for estimating bone mechanical alterations due to disease or drug intervention [1,2]. However, large-scale models require significant computing resources, both in terms of memory and CPU time. Workstations with up to 4GB of RAM and up to 8 processors are becoming widely available, enabling FE simulations without the need for super-computing resources. Here we describe an optimized finite-element solver designed specifically for use with in-vivo MR images of trabecular bone, with the objective of running these simulations within the limitations of standard workstations. The software was applied to compute whole-section axial stiffness in specimens of the human distal tibia from high-resolution MR and microCT images as well as in-vivo 3T MR images in selected subjects being studied as part of a drug intervention trial. Methods: The algorithm is designed on the basis of grayscale bone volume fraction (BVF) maps, where voxels with bone occupancy fraction greater than 20% (0 for pure marrow and 100% for pure bone) were converted into a finite element mesh of hexahedral elements. Each element’s Young’s modulus was set proportional to its BVF assuming a Young’s modulus of 15 GPa and Poisson’s ratio of 0.3 for pure bone. A linear system obtained by minimizing the total elastic strain energy was then solved for the unknown displacements using a preconditioned conjugate gradient (PCG) algorithm. Finally, the bone mechanical parameters were obtained by computing the stress over the strain and a strain-energy map was obtained by calculating the strain energy applied to each voxel. Although the linear system is sparse, the maximum size of the systems being solved is constrained by the available computer memory. Memory requirements can become prohibitive if the coefficient matrix is stored, especially for large-scale models. Therefore, we applied an element-by-element (EBE) approach [3] to the FE simulation algorithm, which avoids assembling and storing the entire global coefficient matrix. Instead, only the BVF and four vertex indices (twenty-four for boundary vertices) associated with each element are stored. The memory usage is hence reduced by a ratio of approximately 3-4. Furthermore, the algorithm takes advantage of the EBE data structure to distribute the collection of elements evenly among different processors, facilitating parallel computing and saving CPU time. The optimized algorithm was applied to different studies including ex-vivo micro-CT images and in-vivo MR images. The in-vivo MR images with 160 μm isotropic voxel size for a 80 x 64 x 10 mm volume of distal tibia of selected subjects were acquired using a 3D FLASE (fast large angle spin echo) pulse sequence [4,5] and a 4channel receive-only phased-array RF coil (Insight, Worcester, MA) on a 3T Siemens Tim Trio scanner (Erlangen, Germany). To obtain the grayscaled BVF maps, each MR image reconstructed from motion-corrected k-space data was first manually masked to isolate the tibial trabecular bone region. Then a local thresholding algorithm [6] was applied to correct for signal variations caused by inhomogeneous sensitivity of the receive coil. The resulting images were then normalized and inverted to obtain the gray-scale BVF maps. To demonstrate the feasibility of analyzing MR images with large number of voxels, the in-vivo MR images were also sinc-interpolated by factors of 3 and 5, yielding 53.3 μm and 32 μm isotropic voxel sizes, respectively. Results and Discussion: Figure 1 illustrates the memory usage (in GB) of the simulations versus number of elements, where a linear relationship was observed. For instance, a desktop computer with up to 4GB RAM can simulate a system with up to about 3.5 million finite elements using the optimized algorithm. Simulation times (in minutes) with 1, 2, 4 and 8 processors, respectively, versus number of elements are given in Figure 2, where all functions fitted are second order polynomials. Parallel computing with 8 processors brought the simulation time down to around 22 minutes (compared with around 2 hours for a single processor) for a system with 12 million elements. Simulations for all non-interpolated in-vivo MR images (160 μm isotropic resolution) contained less than 1 million elements and were completed within 1 minute. In addition to number of processors and number of elements, we noticed that the structure of bone also affected simulation time (i.e., by affecting the number of iterations required). For example, time spent on whole-bone simulation was less than that spent on cuboid subvolume simulation, with both having the same number of elements. The optimized algorithm was also able to perform simulations on interpolated MR images with isotropic voxel sizes of as small as 32 μm, in addition to performing simulations on in-vivo MR images. Figure 3a shows the average intensity projections of the BVF maps of a subvolume of a patient’s distal tibia in coronal view. The average intensity projections of strain-energy (in N•mm) maps of the same subvolume after applying simulated compression to the whole distal tibia along the longitudinal axis are given in Figure 3b (again in coronal view), illustrating that interpolation provides a sharper depiction of the strain map, where each strain trace can be clearly recognized. From comparing the BVF map and the strain-energy map it is also clear that not all trabeculae participated in the load transfer (an example of a region with a significant volume of bone but very low strain is circled in Figure 3). Furthermore, most of the horizontal trabeculae were not loaded because the compression was applied along the vertical direction.